#include "base_ft.h"
#include "lime.h"
#include "read_lime_propagator.h"


int main(int argc,char **argv)
{
   //////////////////////
  // Input
  //////////////////////////


  char kappa_sea[1024];
  int conf;
  int ir;
  int im;

  int  rotate=0;
  int sign=-1;

  char dir_input[1024];
  char dir_output[1024];

  char inputfile_ini[1024];
  char inputfile_end[1024];
  char outputfile_ini[1024];
  char outputfile_end[1024];

	// Default
	sprintf(dir_input,"/home/nucave/workdir/lattice/RENORMALIZATION/DATA/INPUT/TwistedMassNf4/b195/2m/s0-2/datlat/");
        sprintf(dir_output,"/home/nucave/workdir/lattice/RENORMALIZATION/DATA/INPUT/TwistedMassNf4/b195/2m/s0ft-2/");

	conf=501;
	ir=0;
	im=0;
	L=24; T=48;

          ////// Parsing command line

          for(int i=1; i<argc; i++) {
            if(strncmp(argv[i],"-kappa_sea",10)==0)  sscanf(argv[i],"-kappa_sea=%s",kappa_sea);
            if(strncmp(argv[i],"-conf",5)==0)  sscanf(argv[i],"-conf=%d",&conf);
            if(strncmp(argv[i],"-ir",3)==0)  sscanf(argv[i],"-ir=%d",&ir);
	    if(strncmp(argv[i],"-im",3)==0)  sscanf(argv[i],"-im=%d",&im);
            if(strncmp(argv[i],"-dir_input",10)==0)  sscanf(argv[i],"-dir_input=%s",dir_input);
            if(strncmp(argv[i],"-dir_output",11)==0)  sscanf(argv[i],"-dir_output=%s",dir_output);
            if(strncmp(argv[i],"-inputfile_ini",14)==0)  sscanf(argv[i],"-inputfile_ini=%s",inputfile_ini);
            if(strncmp(argv[i],"-inputfile_end",14)==0)  sscanf(argv[i],"-inputfile_end=%s",inputfile_end);
            if(strncmp(argv[i],"-outputfile_ini",15)==0)  sscanf(argv[i],"-outputfile_ini=%s",outputfile_ini);
            if(strncmp(argv[i],"-outputfile_end",15)==0)  sscanf(argv[i],"-outputfile_end=%s",outputfile_end);
            if(strncmp(argv[i],"-T",2)==0)  sscanf(argv[i],"-T=%d",&T);
            if(strncmp(argv[i],"-L",2)==0)  sscanf(argv[i],"-L=%d",&L);

          }

   char filename_input[1024];
   sprintf(filename_input, "%s/%04d/out_r%d_im%02d",dir_input,conf,ir,im);

  char ft_output[1024];
  char sigma_output[1024];
  sprintf(ft_output, "%s/ft_ir%d_im%02d_%04d.out",dir_output,ir,im,conf);
  sprintf(sigma_output, "%s/ff_ir%d_im%02d_%04d",dir_output,ir,im,conf);

          printf("==============================\n");
	  printf("filename_input: %s\n",filename_input);
          printf("ft_output: %s\n",ft_output);
          printf("sigma_output: %s\n",sigma_output);
          printf("L=%d T=%d\n",L,T);
          printf("==============================\n");

  /////////////////////////////////////////
  // end input
  ////////////////////////////////////////


  // momenta
  char momentafile[1024]; sprintf(momentafile,"mom.dat");
  read_momenta_data(momentafile);

  // source position
  int xsource[4];
  char input_source[1024];  sprintf(input_source,"source.dat");  FILE *sfile=fopen(input_source,"r");
  fscanf(sfile,"%d %d %d %d", &xsource[0],&xsource[1], &xsource[2],&xsource[3]);
  fclose(sfile);


  // check endianess

  check_endianess();
  if (big_endian) printf("ATENTION endianness will be transformed \n");

   //////////////////////////////////////////
   // Definitions of the prop and its ft
   ///////////////////////////////////////////
  // PROP IN SPACE 
  colorcolorspinspin *S=malloc(sizeof(colorcolorspinspin)*(T*L*L*L));
  // FOURIER TRANSFORM
  colorcolorspinspin *SFT=malloc(sizeof(colorcolorspinspin)*(Nmom)); 

   // read propagator in position space
   read_binary_colorcolorspinspin(S,filename_input);

  //Some checks
  //print_propagator(S);
  //make_pion(S);

  make_ft_colorcolorspinspin(S,SFT,xsource,sign);

  if (rotate) for (int imom=0; imom<Nmom; imom++) rotate_colorcolorspinspin_to_physical_basis(SFT[imom],!ir);

  //form factors
  complex **s0_ff=malloc(sizeof(complex*)*(Nmom));
  for (int imom=0; imom<Nmom; imom++) s0_ff[imom]=malloc(sizeof(complex)*(3));

  for (int imom=0; imom<Nmom; imom++)
           compute_s0_ff(SFT[imom],s0_ff[imom],SinP[imom],SinP2[imom]);


      // Write ft

      printf("Writing output in %s...\n",dir_output);
      FILE *fout=NULL;
      fout=fopen(ft_output,"wb+");
      if(fout == (FILE *)NULL) {
        fprintf(stderr,"Unable to open file %s for writing\n",ft_output );
        return EXIT_FAILURE;
      }

      for (int imom=0; imom<Nmom; imom++)
        for (int ic1=0; ic1<3; ic1++) for (int is1=0; is1<4; is1++) for (int ic2=0; ic2<3; ic2++) for (int is2=0; is2<4; is2++) 
        {
			double real=creal(SFT[imom][ic2][ic1][is2][is1]);
			double imag=cimag(SFT[imom][ic2][ic1][is2][is1]);
                        fwrite(&real,sizeof(double),1,fout);
                        fwrite(&imag,sizeof(double),1,fout);
         }

     fclose(fout);


      char s0ff_output[1024];
      for (int sigma=0; sigma<3; sigma++)
      {

        // ff S(p)

         sprintf(s0ff_output,"%s_sigma%d.out",sigma_output,sigma+1);
         fout=fopen(s0ff_output,"w");
          if(fout == (FILE *)NULL) {
          fprintf(stderr,"Unable to open file %s for writing\n",ft_output );
          return EXIT_FAILURE;
          }

        for (int imom=0; imom<Nmom; imom++)
		fprintf(fout,"%8.5f %8.10f + I %8.10f\n",SinP2[imom],creal(s0_ff[imom][sigma]), cimag(s0_ff[imom][sigma]));
	fclose(fout);

      }

  return 0;
}
